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We present the first quantitative comparison of two independent general-relativistic hydrodynamics codes, 
the Whisky code and the SACRA code. We compare the output of simulations starting from the same initial 
data and carried out with the configuration (numerical methods, grid setup, resolution, gauges) which for each 
code has been found to give consistent and sufficiently accurate results, in particular in terms of cleanness of 
gravitational waveforms. We focus on the quantities that should be conserved during the evolution (rest mass, 
total mass energy, and total angular momentum) and on the gravitational-wave amplitude and frequency. We find 
that the results produced by the two codes agree at a reasonable level, with variations in the different quantities 
but always at better than about 10%. 
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I. INTRODUCTION 

Given the absence of astrophysically relevant exact solu- 
tions in general relativity and the difficulty to compare results 
from numerical-relativity codes with empirical observations 
(or experiments), it is necessary to find alternative ways to as- 
sess the capacity of existing codes to faithfully describe the 
physical phenomena that they are supposed to simulate, and 
to check the validity of their results. Among the strategies to 
achieve such a reassuring confirmation, the most widely used 
are convergence tests and checks of the violations of the phys- 
ical constraints imposed by the Einstein equations, in partic- 
ular of the so-called Hamiltonian and momentum constraints, 
dictated by the choice of the Arnowitt-Deser-Misner (ADM) 
formalism as a basis for numerical simulations (see Eqs. ( fTOl - 
fTTT i and, e.g., Refs. lll-Ql). Another way to increase the prob- 
ability of having computer codes free from implementation- 
errors and unaffected by possibly wrong and maybe hidden 
assumptions is the comparison of the results of codes inde- 
pendently developed by separate individuals or groups. 

Since 2005, the year of the breakthrough in numerical rel- 
ativity ||3-|al that made it possible to calculate the late inspi- 
ral, merger, and ringdown of a black-hole binary system in 
full general relativity, and to calculate the gravitational waves 
produced in the process, various works (0-121 compared the 
gravitational waveforms computed in vacuum simulations by 
several codes. Their general conclusion is that the available 
codes give consistent results (the difference among codes is 
smaller than the estimated error within each code) and results 
that are good enough for being of use in the quest for the 
detection of gravitational waves through presently operating 
laser interferometers Ill0l - [l2ll or planned detectors Ill3l[l4ll . 

In the present work, with in mind the goals delineated 
above, we perform and publish for the first time a comparison 
between the results of two independent finite-difference codes 
solving the general-relativistic hydrodynamics equations and 
the Einstein equations: the Whisky code lllSl - UTIl and the 
SACRA code yl]. We include in the comparison also impor- 



tant quantities not directly related to the gravitational wave- 
forms, but, while giving a first glimpse of the comparison 
of the wave properties, we postpone to a future article lUSll . 
which may involve a larger number of codes, the detailed anal- 
ysis of the usefulness of the computed waveforms for current 
gravitational-wave detectors. 

We also restrict our attention to the modeling of a single 
physical system, the orbital inspiral of two neutron stars (NSs) 
in irrotational configuration. This system is however one of 
the most promising candidates for early detection of gravita- 
tional radiation and it is seen as the most likely scenario lead- 
ing to the formation of a black hole surrounded by a massive 
torus with properties suitable for being the engine powering 
short-hard gamma-ray bursts Ill9ll . 

We use a spacelike signature (—,+,+,+) and a system of 
units in which c — G — Mq = 1 (unless explicitly shown 
otherwise for convenience). Greek indices are taken to run 
from to 3, Latin indices from 1 to 3, and we adopt the stan- 
dard convention for the summation over repeated indices. 



II. MATHEMATICAL AND NUMERICAL SETUP 

All the details on the mathematical and numerical setup 
used by the two codes have been discussed in depth in pre- 
vious works [2, .3, 20!|. In what follows, we limit ourselves to 
a brief overview, while spelling out the differences between 
the two codes. 

The differences in the implementation of the Einstein and 
hydrodynamics equations between Whisky and SACRA are 
summarized in Table H] 



A. Evolution system for the fields 

We evolve the Einstein equations in the Baumgarte- 
Shapiro-Shibata-Nakamura (BSSN) formalism [21-24]. 

For the Whisky simulations, all the equations discussed 
in this section and in the next are solved using the C CAT IE 



TABLET: Differences between Whisky and SACRA in the schiemes fertile evolution of the spacetime and of the hydrodynamics. See text for 
definitions and further explanations. 





Whisky 


SACRA 


conformal factor 4> 


evolve 4> 


evolve X = e^^** 


primitive matter variables 


p, v\ £ 


p, Ui, £ 


evolved matter variables 


D, S,, T 


D, S„ E 


reconstructed matter variables 


primitive variables: p, u\ £ 


D,u, = Si/D,£ 


local Riemann solver 


Marquina flux formula 


central scheme (Kurganov and Tadmor) 


atmosphere treatment 


constant rest-mass density 


exponentially decreasing rest-mass density 



tion equations for the listed variables. They are: 

[dt - Cp) 7y = -2aA,j , (4) 



code, a three-dimensional finite-differencing code based on 
the Cactus Computational Toolkit lIZSIl . A de- 
tailed presentation of the code and of its convergence prop- 
erties has been presented in 12011 . For tests and details on 
SACRA, see instead iQ]. 

In the BSSN formalism, the spacetime is first decomposed 
into three-dimensional spacelike slices, described by a met- 
ric 7y, an extrinsic curvature A'y, and the gauge functions a (5^ _ £^) ^^^ ^ e'^'^'l-ViVja + a{Rij - SnSi- 



{dt- Cp) (l) = --aK , or {dt - Cp) x= ^axK , (5) 



^TF 



(lapse) and /3^ (shift) (see Sec. lIIBI for details on how we treat 
gauges and ll26ll for a general description of the 3 + 1 split). 
The standard 3h-1 formulation is then modified by introducing 
different variables as follows. The three-metric 7^ is confor- 
mally transformed via 



-^a{KA,j-2A,kA^j), 
{dt - Cp) K = -V'V.a 



A,,A'^ + -K^ + 47r{p,^,,+S) 



(6) 



(7) 
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■ In det 



lij: 



1^J 



6 lij J 



(1) 



and the conformal factor (in Whisky/CCATIE) or a func- 
tion of it (x = e"^"^, in SACRA) is evolved as an independent 
variable, while 7ij is subject to the constraint dct7y = 1. 
The extrinsic curvature is subjected to the same conformal 
transformation and its trace K is evolved as an independent 
variable. That is, in place of Kij we evolve: 



K = tYK,,^-f''K,. 



i,, = e-^^{K,, ~ -j,,K), 



(2) 



with tr Aij = 0. Finally, new evolution variables 



f' = f'T], = -fr 



(3) 



are introduced, defined in terms of the Christoffel symbols of 
the conformal three-metric. 

The Einstein equations specify a well-known set of evolu- 



dtr 



+ |f'aj/3^' - 2A''dja + 2a{rjkA^'' + 6A''dj 



;f^djK-8TTf^Sj), 



(8) 



where Rij is the three-dimensional Ricci tensor, Vi the co- 
variant derivative associated with the three metric 7^, "TF" 
indicates the trace-free part of tensor objects, 5 = ^^^Sij, and 
Padm' '^i' ^iid Sij are the matter source terms defined as 



S, = -j,c.npT''^ , 



(9) 



where ria = (— a, 0,0, 0) is the future-pointing four-vector 
orthonormal to the spacelike hypersurface and T"'^ is the 
stress-energy tensor for a perfect fluid {cf. Eqs. Il27] ll. The 
Einstein equations also lead to a set of physical constraint 
equations that are satisfied within each spacelike slice, 

n = i?(3' +K^ - K,jK'^ - 167rp^oM == 0, (10) 
M' = Vj{K''^ -f^K)-8-KS' = Q, (11) 

which are usually referred to as Hamiltonian and momentum 
constraints, respectively. Here R''^^ = Rijl^-' is the Ricci 
scalar on a three-dimensional timeslice. Our specific choice 



of evolution variables introduces five additional constraints, SACRA the AH is located as reported in ^. 



det7.y = 1, 
trAij = 0, 



(12) 
(13) 
(14) 



Our codes actively enforce the algebraic constraints (fTSl i 
and ( fTSI ). Specifically, after every time evolution, we perform 
a reset as follows: 



Ty ^ [det(7,,)] '/%,, 



(15) 

A,, ^ [det(7,,)]-i/34- - i%-Tr(i,,), (16) 

(17) 



A'^/v +Tr(A„-). 
In SACRA, the additional resetting 



(18) 



is performed. We note that in these adjustments jij and Kij 
are unchanged. 

The remaining constraints, "H, A4^, and (fl4] i. are not ac- 
tively enforced and can be used as monitors of the accuracy 
of our numerical solution. See |27il for a more comprehensive 
discussion of these points. 



B. Gauges 

We specify the gauge in terms of the standard ADM lapse 
function, a, and shift vector, /?' ||28I1 . We evolve the lapse 
according to the "1 + log" slicing condition ll29ll : 



dta - (3'd,a 



-2aK. 



(19) 



The shift is evolved using the hyperbolic F-driver condi- 
tion llll 






I- 



r]B' 



(20) 
(21) 



where ?/ is a parameter which acts as a damping coefficient. 
We set it to be constant and sa 3/A/b, where Afb is the 
baryon mass of one of the stars (for the simulations made with 
Whisky in the present work, the results do not change appre- 
ciably if ?7 is changed at least within a factor 2 of the above 
value). The advection terms on the right-hand sides of these 
equations have been suggested in Il30l - [32ll . 



C. Apparent horizons and gravitational waves 

After the merger, the apparent horizon (AH) formed during 
the simulation is located every few timesteps during the 
evolution. In Whisky this computation is performed both 
with the AHFinderDirect code of 11331 |34[] and in the 
isolated and dynamical-horizon frameworks Il35l - t39ll . In 



For the results reported in the present work, both codes 
extract the gravitational waves using the Newman-Penrose 
formalism, which provides a convenient representation for 
a number of radiation-related quantities as spin-weighted 
scalars. In particular, the curvature scalar 

-Ca/3-ysn"fh'^n'fh^ (22) 



*4 



is defined as a particular component of the Weyl curvature ten- 
sor Cap^s projected onto a given null frame {I, n, m, rh} and 
can be identified with the gravitational radiation if a suitable 
frame is chosen at the extraction radius. In practice, we define 
an orthonormal basis in the three-space {r,0,(p), centered on 
the Cartesian origin and oriented with poles along z. The nor- 
mal to the slice defines a timelike vector t, from which we 
construct the null frame 

l = —{i~r), „= (t + f), m = — (0-i0). 

(23) 
We then calculate ^4 via a reformulation of ( l22b in terms of 
ADM variables on the slice iH: 



^4 = Cij-m^m^ , 



(24) 



where 



Ci. 



Rij — KKij 



K.^Kkj - ie,^'^iK,k (25) 



and eijk is the Levi-Civita symbol. The gravitational-wave 
polarization arnplitudes h+ and /ix are then related to ^1*4 by 
time integrals ll4lll : 



i/i> 



*4 



(26) 



where the double overdot stands for the second-order time 
derivative. Caution should be taken when performing such 
integrals iH]. 

For the extraction of the gravitational-wave signal, both 
codes also implement an independent method, which is based 
on the measurements of the nonspherical gauge-invariant met- 
ric perturbations of a background spacetime II43I1 . The wave 
data obtained in this way give results compatible with the ones 
obtained with the Newman-Penrose formalism and are not re- 
ported here. 



D. Evolution system for the matter 

Both codes adopt a flux-conservative formulation of the hy- 
drodynamics equations ll44l - l46ll . in which the set of conserva- 
tion equations for the stress-energy tensor T^" ~ phu^u'^ + 
pg^'^ and for the matter current density J^ = pu^ (see below 
for definitions), namely 



V^T^'' = 



V^ J'^ = , 



(27) 



is written in a hyperbolic, first-order, flux-conservative form 
of the type 



9tq + 9,f«(q) = s(q) 



(28) 



where f ^'^ (q) and s(q) are the flux vectors and source terms, 
respectively IMTIl . Note that the right-hand side (the source 
terms) does not depend on derivatives of the stress-energy ten- 
sor Furthermore, while the system (|28] | is not strictly hyper- 
bolic, strong hyperbolicity is recovered in a flat spacetime, 
where s(q) = 0. 

The primitive hydrodynamical variables are the rest-mass 
density p, the specific internal energy e measured in the rest- 
frame of the fluid, and the fluid three-velocity (defined as 
v'^ = u^/W + /3'^/a (contravariant components) in Whisky 
and as u^ (covariant components) in SACRA, where u^ is the 
four-velocity measured by a local zero-angular-momentum 
observer; SACRA defines contravariant components of the 
three-velocity as V^ ~ u^ /u°). The Lorentz factor is defined 
as 






(29) 



There is then an equation of state (EoS) relating pressure, rest- 
mass density and internal-energy density. 

Following 114511 . in order to write system ( |27] | in the form 
of system ( l28T l, the primitive variables are mapped to a set of 
conserved variables q = (D, Si, E) via the relations 

D = ^W p = e^^'I'W p , 



Si = Dui = y/jphW'^Vi 



(30) 



E 



/7 {phW"^ 



P 



D = Di, 



where h = l + e+p/pis the specific enthalpy, Ui = hui is the 
specific momentum, and e = hW — p/{pW) is the specific 
energy. 

In this approach, all variables q are represented on the nu- 
merical grid by cell-integral averages. The functions that the 
variables q represent are then reconstructed within each cell, 
usually by piecewise polynomials, in a way that preserves 
conservation of the variables q ll48ll . This operation pro- 
duces two values at each cell boundary, which are then used 
as initial data for the local Riemann problems, whose (ap- 
proximate) solution gives the fluxes through the cell bound- 
aries. A method-of-lines approach ll48ll . which reduces the 
partial differential equations ( |28] | to a set of ordinary differ- 
ential equations that can be evolved using standard numer- 
ical methods, such as Runge-Kutta or the iterative Cranck- 
Nicholson schemes 11491 15011 . is used to update the equations 
in time (see lUSll for further details). Here, we employ the 
4*^-order Runge-Kutta method (see below). 

Various reconstruction methods are implemented in 
Whisky and SACRA, but here we always use the piecewise 
parabolic method (PPM) Iplll . Both codes implement the 
scheme of Kurganov-Tadmor 115211 (which is a variation of 
the HLLE approximate Riemann solver ll53ll ). but Whisky 
gets better results employing the Marquina flux formula ll54ll 
(see 1 15, 16] for a more detailed discussion). A comparison 
among different numerical methods in binary-evolution 
simulations was reported in 13, QM ■ 

There are differences between Whisky and SACRA in sev- 
eral implementation choices. In Whisky, the variables whose 



evolution is computed are D, Si, and t = E — D. SACRA 
adopts as evolution variables D, Si, and E. Furthermore, the 
PPM reconstruction is performed by SACRA on the variables 
p, Ui ~ Si/D, and e, while Whisky reconstructs the primi- 
tive variables p, w\ and e. 

Other differences are present in the conversion from the 
evolved conservative variables back to the primitive variables, 
which are used to calculate the fluxes and the source terms of 
the equations. Such a conversion cannot be given in an ana- 
lytical closed form (except in certain special circumstances). 

Whisky implements the following procedure to do the 
conversion. One writes an equation for the pressure 



P^ 



-p[p(q,p),£(q,p)] =0 



(31) 



where p is the value of the pressure to be found and 
p[p(q, p), e(q, p)] is the pressure as obtained through the EoS 
in terms of the updated conserved variables q and of p itself. 
This is done by inverting ( [30l ) to express p and e in terms of 
the conserved variables and of the pressure only: 



D 



P = 



P + D 



^iT+p + D)^-S^, 



e = D- 



^(T+p + Df- S'^ -pW -D 



where 



W 



D 



V^ 



P 



D) 



52 



(32) 



(33) 



(34) 



is the Lorentz factor, expressed in terms of the conserved vari- 
ables, and 



52 = j'^S^S, 



(35) 



Then (l3Tl i is solved numerically. In Whisky we use a 
Newton-Raphson root finder, for which we need the deriva- 
tive of the function with respect to the dependent variable, i.e. 
the pressure. This is given by 



-I 
dpi 



P 



p[/o(q,p),£(q,p)]} 

dp{p, e) dp dp{p, e) de 



1 



where 



dp 

dp 

de 
dp 



dp dp 



DS^ 



de dp ' 



V^ 



Df-S'^{T- 
pS^ 



D) 



p[{t+p + D)^-S^]{t+p + D)'' 



(36) 

(37) 
(38) 



and where dp/ dp and dp/de are given by the EoS. Once the 
pressure is found, the other variables follow simply. 

In SACRA, the conversion is performed in the following 
way. From the normalization relation u'^u^ = —l,Wis ex- 
pressed in terms of h and of the evolved values of Mi (= Si/D) 
and 7*-' : 



W^ 



1 



1^'UiUj 
/l2 ■ 



(39) 



For the EoSs chosen in the present work, h is regarded as 
a function of W for the evolved values of D, e, and cj) be- 
cause of the relation e = hW — pe^'^/D and of the fact that 
p is written as a function of h, p{= De'^"^ /W), and W as 
p = p{h, p) = p{h, W). Substituting the resulting relation 
for h = h{W) into Eq. ( |39] |, we obtain a one-dimensional 
algebraic equation for W both for the F-law EoS and the 
piecewise-polytropic EoS (see Sec. HID 2] ). We solve this de- 
rived equation using the Newton-Raphson method, for which 
we need to take a derivative of the equation F{W) = with 
respect to W. This is rather straightforward, and straightfor- 
ward is also the determination of the variables h, e, and P, 
once the equation for W is solved. 



1. Treatment of the atmosphere 

At least mathematically, the region outside our initial stel- 
lar models is assumed to be perfect vacuum. Independently 
of whether this represents a physically realistic description of 
a compact star, the vacuum represents a singular limit of any 
conservative scheme for hydrodynamical evolution and must 
be treated artificially. Both codes follow a standard approach 
in computational fluid-dynamics, that is the addition of a ten- 
uous "atmosphere" filling the computational domain outside 
the star 

Of course, the density of the atmosphere should be as small 
as possible, in order to avoid spurious effects. The evolution 
of the hydrodynamical equations in grid zones where the at- 
mosphere is present is the same as the one used in the bulk 
of the flow. When the rest-mass density in a grid zone falls 
below the threshold set for the atmosphere, that grid zone is 
not updated in time and the values of its rest-mass density, 
internal-energy density, and velocity are set to those of the 
atmosphere. 

Both codes treat the atmosphere as a zero-coordinate- 
velocity perfect fluid governed by a polytropic EoS with the 
same adiabatic index used for the bulk matter, or, in case of the 
piecewise-polytropic EoS (see Sec. lIID2b . the same adiabatic 
index as the one used in the outer parts of the star' . However, 
the values of the rest-mass density assigned to the atmosphere 
are different. In Whisky, the rest-mass density is set to be 
constant and several (10 in the present simulations) orders of 
magnitude smaller than the initial maximum rest-mass density 

Pmax ass. 

In SACRA the rest-mass density is assigned as 



P 



P.. e^-''/''" r>r„. 



(40) 



where p^^^^ — p^^^ x 10^^ is chosen, r^ is a coordinate 
radius of about 1Q-2QM f^^^^, where M^^^^^ is the ADM mass 
of the system. In both codes, also the internal-energy density 
e is then recomputed from p according to the polytropic EoS. 



' In this case, also the polytropic constant K for the atmosphere is chosen to 
be the same as the one in the outer parts of the star. 



For both codes, with such a choice of parameters, the rest 
mass of the atmosphere is at least a factor 10~^ smaller than 
the rest mass of the NSs. Thus, spurious effects due to the 
presence of the atmosphere, such as accretion of the atmo- 
sphere onto the NSs and the black hole, the resulting dragging 
effect against orbital motion, gravitational effects, and effects 
on the formation and dynamics of the disk around the merged 
object play a negligible role in the present context. 



2. Equations of state 

In this work we present results obtained with two EoSs: 
a simple "F-law" or "ideal-fluid" EoS and a piecewise- 
polytropic EoS Ipql . For the ideal- fluid EoS, the pressure is 
given as 



p=(F-l)p£ 



(41) 



where F is the adiabatic index. When using the ideal-fluid 
EoS ( |4TI ). nonisentropic changes can take place in the fluid 
and, in particular, shocks (which are always present in the 
mergers and which may play important roles) are allowed to 
transfer kinetic energy to internal energy. On the other hand, a 
carefully chosen piecewise-polytropic EoS may mimic more 
closely a realistic EoS. The parametrised EoS we consider 
consists of two polytropes interfacing at a density pQ. The 
relations between the hydrodynamical quantities are {i = 
0,1)111 



P 



K,p' 



K, 



(1 + «.)P+JT^/- 



(42) 
(43) 



where Ki are the polytropic coefficients, and Fj are the poly- 
tropic exponents in the different intervals of rest-mass density. 
Furthermore, the constants ai, which guarantee continuity, are 



ao = , 

g(Po) 
Po 



ai 



- 1- 



Po 



Ti-l' 



(44) 
(45) 



In our simulations we used the parameters of model B of IS 
namely: 



Po = 1.630497500125504 x lO""* g cm 



(46) 



< p < Po : To 

A'o 

p> Po ■ Ti 

ai 



1.35692395 , 

0.35938266 x lO^'' cgs units, 

3.0, 

0.15982116 X 10"^ cgs units, 

0.01088158737430845 . 



In the presence of shock heating, part of the kinetic en- 
ergy is converted into thermal energy. To model this property, 
the original piecewise-polytropic EoS is modified by adding a 
thermal contribution to the pressure 



Pf 



th 



(Fth- l)p(e-eo). 



(47) 



TABLE II: Differences in the implementation of tlie AMR of SACRA and Whisky. See text for definitions and furtlier explanations. 





Whisky 


SACRA 


prolonged and restricted variables 


conserved variables: D, Si,T 


D,u,{^S,/D),h 


interpolation for the prolongation 
of the hydrodynamical variables 


ENO: 3"^ order in space, 2'"* order in time 


Lagrangian: S**" order in space (reduced to 1°' order in 
case of failure), 2"'^ order in time 
(reduced to 1°* order at extrema) 


buffer zones 


12 


6 


overlapping same-level grids are 


evolved as a single grid 


evolved independently (but using the average of the values 
of the two grids at overlapping points) 



where Fth is the adiabatic index for this correction and Eq is 
given by Eq. ( |43] l. In the absence of shocks, e is equal to Eq 
and thus Pth = 0. In the simulations of this work SACRA has 
adopted Fth = Fq while in Whisky the thermal correction 
was not applied (Pth = ). As the figures of this work show, 
at least in the inspiral phase the difference in the adopted EoS 
does not have an influence. 



E. Adaptive Mesh Refinement 

There are similarities and differences in the implementation 
of the adaptive mesh refinement (AMR) in the two codes. In 
the following we spell them out in detail. 

Both codes employ a vertex-centered Berger-Oliger 115811 
mesh-refinement scheme adopting nested grids with a 2 : 1 re- 
finement factor for successive grid levels. In the simulations 
made for the present work both codes used a set of coarser 
fixed grids and finer moving grids, centered around each 
star Whisky makes use of the Carpet mesh-refinement 
driver 1521 ■ The higher-resolution moving grids are centered 
around the local maximum in the rest-mass density p of each 
star In SACRA, instead, the grids are centered around the lo- 
cal maximum of the conserved variable D. 

Both codes employ centered 4"^-order finite-differencing in 
space for evaluating spatial derivatives of the geometric quan- 
tities, except for the shift advection terms that are calculated 
with upwinding derivatives to improve accuracy. For the time 
integration, the 4*'^-order Runge-Kutta scheme is adopted. To 
evolve quantities near the refinement boundaries of a refined 
grid, both codes introduce buffer zones, where the variables 
are computed ("prolonged" and "restricted") in a special way 
and not with the time-update scheme used for all other non- 
refinement-boundary points. 

In Whisky/Carpet, we use 12 buffer points, 3 for each 
substep of the adopted time-integration scheme. The values 
of the needed quantities at the buffer points are computed 
from the coarser grid through interpolation as follows: For 
the spacetime variables, 5"^-order Lagrangian interpolation 
in space and 2"'^-order Lagrangian interpolation in time are 
used; For the hydrodynamical variables, 3"^-order ENO ll60ll 
interpolation in space and 2"'^-order ENO interpolation in 
time are used. The prolonged and restricted variables are the 
conserved evolved ones: D, Si, and r. The interpolation is 
done whenever the first Runge-Kutta time integration is being 
carried out. 



In SACRA both prolongation and restriction are carried out 
on D, Ui{~ Si/D), and h. Following ll6lll . 6 buffer points 
are introduced. The quantities at the buffer zones are pro- 
vided from the corresponding coarser domain by the follow- 
ing procedure. For space interpolation, 5"^-order centered La- 
grangian interpolation in space is carried out using the nearby 
6 points of the coarser grid. This is done both for spacetime 
and hydrodynamics variables. For the latter, this interpolation 
scheme could fail, in particular in the vicinity of the surface of 
the stars, where D is small and varies steeply. The reason for 
this possible failure is that the interpolation may give a nega- 
tive, and so unphysical, value of D or h — 1. If the 5"^-order 
Lagrange interpolation produces D < Z^min or h < 1, F^'- 
order (i.e., linear) interpolation is adopted. Z^min is chosen 
to be Dmax/lO^, where I?inax is the initial maximum value 
of D. Linear interpolation cannot be used in general for all 
points because it is too dissipative. As in Whisky the in- 
terpolation is also done whenever the first Runge-Kutta time 
integration is being carried out. 

For the update of the buffer zones SACRA implements, in- 
stead, the following: i) For the inner three buffer points all the 
quantities are evolved using the 4"^-order finite-differencing 
scheme. Since there is a sufficient number of buffer points to 
solve the evolution equations in the inner three buffer points, 
no interpolation is necessary; ii) For the fourth buffer point, all 
the quantities are evolved using a 4"^-order finite-differencing 
scheme with no interpolation, except for the transport terms 
for the geometry such as /S'^d^jij, for which 2"'^-order finite 
differencing is employed when /?'"' has an unfavorable sign; 
iii) For the two outer buffer points, 2"'^-order Lagrangian in- 
terpolation in time of the coarser-grid quantities is carried out. 
This time-integration procedure is applied to both spacetime 
and hydrodynamical variables, but for the latter there is an 
additional check. 

The interpolated value at a finer-grid time step is obtained 
from the values at the three time levels of the coarser grid, 
say, n — 1, ?7, and n + 1 (note that n does not denote the 
Runge-Kutta time step). The interpolation is necessary for 
determining the values at a time t that satisfies t" < t < <"+^. 
Defining Q as D, Ui, or h, and Q" as the value of the variable 
Q at time t", SACRA checks whether (Q"+i - Q'^){Q"- - 
Qn-i-j ^ Q ^jjj j£gQ adopts P* -order interpolation, using only 
Qn+i ^jjj Qn Namely, a limiter procedure is introduced. 
This robust prescription provides numerical stability yy. 

The two domains in the finer levels often overlap. In such 
cases, the values of all quantities should agree with each other. 



TABLE III: Properties of the initial data: proper separation between the centers of the stars d/A/^^^j ; baryon mass Mb of each star in units 
of solar mass; total ADM mass A/^^^, in units of solar mass, as measured on the finite-difference grid; total ADM mass M^^^-^ in units of 
solar mass, as provided by the Meudon initial data; angular momentum J, as measured on the finite-difference grid; angular momentum J, 
as provided by the Meudon initial data; initial orbital angular velocity Jlo; mean coordinate equatorial radius of each star r^ along the line 
connecting the two stars; maximum rest-mass density of a star Pmax- The columns for Maom and J contain the value for Whisky (left) and 
the one for SACRA (right). Note that the values of A/adm and J are computed through a volume integral in Whisky, while in SACRA they 
are computed through the extrapolation to r —i- cxd of the ADM masses and angular momenta calculated as surface integrals at finite radii r. 



EoS for the model 


rf/^AOM 


Mt 
(Afs) 


(Ms) 


Maom 

(Mo) 


J 

(xlO'^'gcmVs) 


J 

(xlO*»gcmVs) 


(rad/ms) 


(km) 


Pmax 

(g/cm^) 


Ideal fluid (T = 2) 
Piecewise polytropic 


12.6 
15.4 


1.779 
1.502 


3.251,3.256 
2.676, 2.680 


3.233 
2.668 


8.921,8.930 
6.492, 6.506 


8.922 
6.491 


1.906 
1.664 


12.23 
8.48 


7.58 X 10" 
9.77 X 10" 



but, since in SACRA the evolution equations for the two do- 
mains are solved independently, the values do not always 
agree exactly. Let us denote with Qi and Q2 the values on 
the two domains of an individual refinement level. In order 
to guarantee that they are the same, in SACRA the average of 
the two values is used: Qi = Q2 ^ (Qi + (32)/2. When a 
buffer point of one of the two domains overlaps with a point in 
the main region of the other domain, the values at the point of 
the main region are copied to those at the buffer point. When 
two buffer zones overlap at some points, the simple averaging 
described above is again used. 

In Whisky, when domains of the same refinement level 
would overlap, the whole level is automatically resplit in 
(smaller and more numerous) nonoverlapping domains, so in 
practice they continue to be evolved as a single grid, without 
requiring averaging. For more details on the Carpet code 
seelUl. 

For both codes, at the outer boundaries of the coarsest re- 
finement level, an outgoing boundary condition is imposed for 
all the geometric variables. The outgoing boundary condition 
is the same as that suggested by Shibata and Nakamura 1I22II . 
Flat boundary conditions are applied to the matter variables. 

Both codes can add artificial dissipation to the source terms 
of the Einstein equations. In particular, for the schemes pre- 
sented in this work, they could use 5"^-order Kreiss-Oliger- 
type dissipation ||62I1 as Q; — > Q; — crhfQl where Qi is a 
quantity in the ^th level, hi is the spacing of the Z— th level, 
Ql is the sum of the sixth derivatives along the x, y, and z 
axis directions, and a is a constant of order 0.1. The results of 
the present work were obtained without artificial dissipation 
for SACRA and with artificial dissipation for Whisky. 

Standard SACRA simulations for NS-NS binaries are per- 
formed with 7 or 8 refinement levels, in particular 3 or 4 
coarser levels composed of one domain and 4 finer levels com- 
posed of two domains. The time step for each refinement 
level, dti, is determined as follows: 



dti 



h2/2 ioi-0<l<2 
hi/2 ioi-2<l <L-1. 



(48) 



Namely, the Courant number (expressed in terms of the speed 
of light) is 1/2 for the finer refinement levels with / > 2, 
whereas for the coarser levels, it is smaller than 1/2. The rea- 
son why a smaller Courant number is chosen for the coarser 
levels is that with a Courant number as high as 1/2, numerical 



instabilities occur near the outer boundary. This is an inher- 
ent problem of the adopted F-driver gauge condition [? ] and 
it does not appear in the Whisky simulations of the present 
work only because the resolution in the coarsest grids is still 
high enough. 

In standard Whisky simulations for binary systems 6 re- 
finement levels are used, the two finest of which move follow- 
ing the stars. In addition to the moving grids, a set of refined 
but fixed grids is set up at the center of the computational do- 
main so as to capture the details of the Kelvin-Helmholtz in- 
stability (see 121]). The Courant factor is 0.35 for all levels. 

In the Whisky simulations for the present work, a re- 
flection symmetry condition across the 2; = plane and a 
TT-symmetry condition^ across the x = plane are used, 
while SACRA adopts only the reflection symmetry across the 
z plane. 

The differences in the implementation of AMR between 
SACRA and Whisky are summarized in Table HIl 



F. Initial data 

The initial configurations for our relativistic-star binary 
simulations are produced using the multidomain spectral- 
method code, LORENE, which was originally written by the 
group working at the Observatoire de Paris -Meudon 11631 [6411 
and which is publicly available ll65ll . Specific routines are 
used to transform the solution from spherical coordinates to 
a Cartesian grid of the desired dimensions and shape. 

These initial data, which we refer to also as the "Meudon 
data", are obtained under the assumptions of quasiequilib- 
rium and of conformally-flat spatial metric. The initial data 
used in the simulations shown here were produced with the 
additional assumption of irrotationality of the fluid flow, i.e. 
the condition in which the spins of the stars and the orbital 
motion are not locked; instead, they are defined so as to have 
vanishing vorticity. Initial data obtained with the alternative 
assumption of rigid rotation were not used because, differ- 
ently from what happens for binaries consisting of ordinary 



^ Stated differently, we evolve only the region {x > 0, z > 0} applying a 
180° -rotational-symmetry boundary condition across the plane at a; = 0. 



TABLE IV: Properties of the initial grids: number of refinement levels (including the coarsest grid); number of finer levels that are moved to 
follow the stars; spacing of the finest level; length of the side of the finest level; spacing of the coarsest level; outer-boundary location. All 
lengths are expressed in km. HR, MR, and LR denote the high, medium, and low resolutions, respectively. For Whisky the two resolutions 
are in a ratio of 5/4, while for SACRA the ratio between LR and MR is 50/42 and the ratio between MR and HR is 1.16. 



Model 


n° of levels 


n° of moving levels 


finest spacing 


extent of finest grid 


coarsest spacing 


outer-boundary location 


Whisky ideal fluid 


6 


2 


0.1773 


44.33 


5.67 


380 


Whisky piecewise, HR 


6 


2 


0.1773 


44.33 


5.67 


760 


Whisky piecewise, LR 


6 


2 


0.2216 


44.33 


7.09 


760 


SACRA ideal fluid 


8 


4 


0.1387 


6.656 


17.75 


852 


SACRA piecewise, HR 


7 


4 


0.1746 


9.428 


11.17 


603 


SACRA piecewise, MR 


7 


4 


0.2025 


10.13 


12.96 


648 


SACRA piecewise, LR 


7 


4 


0.2411 


10.13 


12.96 


648 



stars, relativistic-star binaries are not thought to achieve syn- 
chronization (or corotation) in the timescale of the coales- 
cence ll66ll . 

The initial models for the binaries have been chosen so as 
to allow significant possibilities of comparison between the 
codes and at the same time to limit the required computa- 
tional time. In particular, after performing several orbits and 
merging, prompt collapse to a Kerr black hole occurs. As said 
above, we chose two EoSs, the ideal-fluid EoS ( 1411 1'^ and a 
piecewise-polytropic EoS ( |43] |. For the latter, the initial data 
have been kindly provided by K. Taniguchi. Note that the 
model with the ideal-fluid EoS has been often used in previ- 
ous work (e.g. lHll). 

Some of the physical quantities of the initial configura- 
tions are reported in Table |III1 In brief, they are equal-mass 
configurations with an initial proper distance between stellar 
centers of about 60 km (initial orbital frequency 0.303 kHz 
and 0.265 kHz, respectively for the ideal-fluid model and for 
the piecewise-polytropic model). The chosen rest masses of 
A/^P = 1.779Mq and M^^ = 1.5O2Af0, respectively for 
the two models, lead - as desired - to prompt collapse to black 
hole. 



G. Specific grid setup for the reported simulations 

For the higher-resolution run with Whisky, the spacing of 
the finest of the six grid levels is h = 0.120 Mq ~ 0.1773 km 
and the spacing in the wave zone (the coarsest grid) is h = 
3.84 M© ~ 5.67km. For the lower-resolution run the spac- 
ing is h = 0.150 M0 ~ 0.2216 km on the finest grid and 
h = 4.80 Mq ~ 7.09 km on the coarsest grid. The finest 
grid always covers the whole stars. For the simulations with 
the ideal-fluid model the outer boundary is located at about 
380km while in the case of the piecewise-polytropic model, 
for both resolutions, the outer boundary is at about 760 km. 
Except for the outer boundary location and the grid spacing, 
the AMR grid structure was the same for all the runs. 



^ The initial data for the simulations adopting the ideal-fluid EoS are set up 
as a simple polytropic EoS with polytropic constant K = 123.6 (in units 
of c = G = Mq = 1). 



For the runs with SACRA, for the ideal-fluid model, the grid 
structure is essentially the same as in t3s]; the finest of the eight 
grid levels has h = 0.0938 Mq ~ 0. 1387 km. For the simula- 
tions with the piecewise polytrope, the computational domain 
is composed of seven grid levels with the finest grid resolu- 
tion being h = 0.1182 Mq ~ 0.1746 km at the high resolu- 
tion, h = 0.1370 M0 ~ 0.2025 km at the medium resolution, 
and h ~ 0.1631 Mq ~ 0.2411km at the lower resolution. 
The resolution in the wave zone (for the coarsest grid level) is 
h ~ 11.17km for the high-resolution run and h ~ 12.96km 
for the others. The boundary of the finest grid is at 60% of 
the stellar radius (along a coordinate axis, att — 0) while the 
second finest grid covers all the stars for the run with the ideal- 
fluid EoS, whereas for the run with the piecewise-polytropic 
EoS, the finest grid covers the stellar radius completely (the 
boundary of the finest grid is at 115% of the stellar radius). 
The outer boundary is at about 852 km for the simulations 
performed with the ideal-fluid EoS and at about603 km or 
648 km for those performed with the piecewise-polytropic 
EoS, for the high-resolution run and the other runs respec- 
tively. 

As already noted in Sec. Ill El another difference between 
the grid setups of the two codes is the adopted symmetry. Both 
codes compute only the z > portion of the {x, y, z} Carte- 
sian coordinate numerical domain, but, while SACRA calcu- 
lates all the z > portion. Whisky calculates only the .t > 
part of the remaining domain, taking advantage of the 180° 
degree rotation symmetry characterising equal-mass binaries. 

The properties of the grids adopted in the simulations with 
the two codes are summarised in Table II VI 

For the setup of the piecewise-polytrope high-resolution 
run. Whisky, which heavily exploits large parallel facili- 
ties, uses approximately 22 x 10® grid points and the to- 
tal required memory for the high-resolution run is about 640 
GBytes. SACRA, instead, which has been specifically devel- 
oped for being able to perform production simulations even on 
a laptop computer, uses about 7 x 10® grid points and about 
11.6 GBytes of memory. For Whisky, the total CPU time 
for the high-resolution piecewise-polytrope run was about 450 
CPU hours on 320 processors of the Ranger cluster (at the 
Texas Advanced Computing Center; the processors are AMD 
Opteron Quad-Core 64-bit, with clock frequency 2.3 GHz) 
and for SACRA it was about 2000 CPU hours on a Quad-Core 
machine of Core-i7X processors with clock frequency 3.33 
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GHz. 



III. COMPARISON OF THE RESULTS 

As also described in ||j,|3|]i the chosen initial data for the 
present study are such that the stars orbit about 3 times and 
7 times, respectively for the two models with different EoSs, 
before merging. As can be seen in Fig. [1] the rest-mass 
density at the stellar centres remains approximately constant 
for the first 6 ms (in the case of the simulation with the 
ideal-fluid EoS) or 15 ms (for the piecewise-polytropic case) 
and then decreases, indicating an expansion of the stars due 
to the tidal force, just before the merger. As expected from 
the (high) mass of the chosen models, the merged object then 
immediately collapses to a black hole and the AH is measured 
for the first time at about 8 and 18 milliseconds, respectively 
for the two models with different EoSs {cf. the highest 
resolutions). The mass 7\/bh and the angular-momentum 
parameter a = ^bh/(^^bh)^ of the resulting black hole 
are measured by both codes. The values after the ringdown 
for the piecewise-polytropic EoS are A/g^^ ^ ~ 2.633A/0, 
^jSACRA ^ 2.637Mo (a relative difference of 0.15%), and 
^Whisky ^ 0.79, a^*™* = 0.80 (a relative difference of 
1.2%). For the ideal-fluid EoS the values of the black hole 
are M""'"^ = 3.22Afm, Ml^^'' = 3.2IM0, and a = 0.84 



'BH 

for both codes 



'BH 



Having briefly summarised the dynamics of the system, 
we present now first a comparison between some quanti- 
ties produced in evolutions performed with SACRA and with 
Whisky, each in what is thought to be a good configura- 
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FIG. I: (Colour online) Comparison of the time evolution of the 
maximum of the rest-mass density for the two models (with different 
EoSs) described in Sec. lIID2l For ease of interpretation, we remind 
the reader that in our units p — 1 x 10""^ corresponds approximately 
to 6.18 X 10^-*g/cm^ 



tion in terms of accuracy, violation of the ADM constraints, 
and cleanness of gravitational waves. Furthermore, for the 
piecewise-polytropic EoS we present for each code results ob- 
tained at two or three resolutions. 

From Fig. [1] one can see immediately that the time of the 
merger depends considerably on the grid resolution, for both 
codes, but in a stronger fashion for Whisky. As is well 
known, the conservation of the angular momentum in nu- 
merical simulations of binary compact objects is a delicate 
issue, which can have very visible effects like the ones in 
Fig. [T] Even if the merger and post-merger dynamics may 
not be sensible to the exact timing of the inspiral, the phase 
of gravitational waves is affected and so this effect must be 
carefully taken into account when producing templates for 
gravitational-wave data analysis. For example, |3] attempted 
to do so by estimating, given a specific initial-data configu- 
ration, the 'real' merger time at infinite resolution through an 
extrapolation based on the results of simulations of the same 
model at different resolutions. Anyway, we are here interested 
in the comparison of the codes and note that, when the dif- 
ferences due to the resolution are subtracted by time-shifting 
the curves, the evolutions of the rest-mass density in the two 
codes are very similar. As said, a proper analysis of the phase 
difference of the gravitational waves from the various codes 
and resolutions will be reported in a future work jlSll . 

We continue the discussion of the results in a more quan- 
titative manner by comparing the time evolution of the rest 
mass, which should be a conserved quantity as no matter is 
seen leaving the numerical domain through the outer bound- 
ary during the simulation. One can see in the left panel of 
Fig.|2]that both codes conserve the rest mass at very high ac- 
curacy, but in SACRA the violation is of the order of 10^'^ 
while in Whisky it is of the order of 10^^. More in detail, 
the dot-dashed black line refers to the high-resolution SACRA 
run, which of course shows an improvement over the medium 
(continuous red line) and lower-resolution ones (dotted green 
line). The convergence is achieved approximately at second 
order. The curves referring to Whisky look constant on the 
main panel, but in the subpanel one can notice the minute 
increase in the rest mass even in the high-resolution results 
(short-dashed blue line). The curve referring to the low reso- 
lution (long-dashed magenta) drops at the time of AH forma- 
tion because the matter inside the horizon is not included in 
the computation of the rest mass. 

The reason of the relatively worse conservation in SACRA 
(as said, the conservation is very good in absolute terms also 
for SACRA) is to be found in the presence of a refinement 
boundary very close to the stellar surface. In the orbital phase, 
oscillations due to the tidal deformation of the NS cause the 
matter to cross the finest refinement level and the small errors 
due to the interpolation in the buffer zones are larger where the 
density is larger Also in Whisky, if a refinement boundary 
is placed inside the stars, the violation of the conservation of 
the rest mass is larger (^ 10""'). 

The right panel of Fig. |2] which - as the left one - refers 
to the piecewise-polytropic EoS, shows then the conservation 
of the energy, namely the sum of the ADM mass computed 
on the numerical domain and of the energy carried by gravita- 
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FIG. 2: (Colour online) Left: Comparison of the time evolution of the rest mass (normalised to the initial value). The inset is a magnification 
of the higher-resolution Whi sky curve, in the form M/Mq — 1. These data refer to the piecewise-polytropic EoS. As explained in the text, the 
larger variations in the SACRA data are due to the choice of grid structure. Right: Comparison of the time evolution of the sum (normalised 
to the initial value) of the ADM mass measured on the numerical grid and the energy carried away from the grid by gravitational waves. This 
quantity should be conserved. These data refer to the piecewise-polytropic EoS. Note that the data for the low resolution of Whisky are not 
reliable after the formation of the AH (t ~ 13.4 ms for this simulation), because the volume integral with which the ADM mass is computed 
contained also the points inside the horizon. See text for more details. 



tional waves outside the numerical domain. Such a quantity, 
normalised to its initial value (the initial ADM mass) should 
be constant and the figure shows the deviation of the results 
from constancy. The colours and line types are the same as in 
the left panel. At the highest resolutions, both Whisky and 
SACRA conserve this quantity very well, at the order of 1 per 
1000 during the inspiral and at better than 1% overall. 

Some of the differences in the curves referring to the two 
codes (in particular the 'smoothness') are due to the differ- 
ent way of computing the ADM mass. Whisky performs a 
volume integral with the formula 



M 



ADM, Vol 



ft 



a\/ll^ 1 {dklji - dajk) 



(fin 



(49) 



and in the simulations of this work it does not exclude the 
points inside the AH from the computation, so the values of 
the ADM mass given by Whisky after the appearance of the 
AH are affected by gross errors. SACRA, instead, uses a sur- 
face integral on a spherical surface far from the central ob- 
jects. This method gives consistent results after the formation 
of the AH, but is more sensitive to small metric oscillations 
in the vicinity of the chosen surfaces, which lie in the coarse 
resolution region; the 'roughness' of the curves follows from 
this. 

Figure [5] which refers to the piecewise-polytropic EoS, 
shows then the conservation of the angular momentum, de- 
fined here as the sum of the angular momentum computed on 
the numerical domain and of the angular momentum carried 
by gravitational waves outside the numerical domain. Such 
a quantity, normalized to its initial value (the initial angular 



momentum) should be constant and the figure shows the devi- 
ation of the results from constancy. The colours and line types 
are the same as in the previous figures. At the highest resolu- 
tion. Whisky conserves this quantity very well, at better than 
1%, and also for SACRA deviations from constancy are of the 
same order, even if larger oscillations are visible. The differ- 
ence in the computation of the angular momentum, analogous 
the one for the ADM mass, is also here at the origin of the 
difference in the smoothness of the curves. Namely, Whisky 
performs a volume integral with the formula 116911 
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and excludes from the integral the points inside the AH. How- 
ever, if the angular momentum of the black hole is added to the 
one computed above, the coiTect time evolution of the quantity 
in Fig.[3]is recovered, except for an interval just after the AH 
formation, when the AH is small and covers only a few grid 
points, and so the measurement of its angular momentum is 
inaccurate. SACRA, instead, uses also here a surface integral. 
As previously noted, also from the time evolutions in Fig. [3] 
one sees that the conservation of the angular momentum at 
these resolutions depends in a stronger way on resolution for 
the Whisky code with respect to SACRA. In addition, one 
can see that, while also the SACRA data show convergence 
almost everywhere, in some time intervals the behaviour at 
different resolutions is not convergent, for example at the 
spike around 2.5 ms. The reason is not completely clear at 
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FIG. 3: (Colour online) Comparison of the time evolution of the an- 
gular momentum, computed as the sum of the angular momentum 
measured on the numerical grid and the angular momentum carried 
away from the grid by gravitational waves. These data refer to the 
piecewise-polytropic EoS. As already noted for Fig. [2] also here the 
data for the low resolution of Whi sky are not reliable after the for- 
mation of the AH (t ~ 13.4 ms for this simulation), because the 
volume integral with which the angular momentum is computed ex- 
cludes the contribution of the black hole. See text for more details. 



the moment, but we think that this is probably related to 
the low resolution of the coarsest grid, where the surface on 
which the angular momentum is computed is located [note 
that accurate extraction of angular momentum requires an 
accurate computation of parts of the extrinsic curvature that 
are 0{r~^) and these are much smaller than the leading-order 
wave part of 0(r~^)]. If the angular momentum is computed 
on surfaces that lie on the finer levels, the differences in the 
wrong direction caused by resolution are much smaller (but 
the value of the angular momentum is less accurate). 



We now proceed to analyze gravitational waves extracted 
from the simulations. The data presented here are extracted 
from the numerical simulations at distances from the origin of 
the axes in the interval 300 ^ 600 km. For building templates 
to be used in the analysis of the data taken by the gravitational- 
wave detectors, the accurate knowledge of the frequency of 
the waves is of special importance. Thus we first show in 
the left panel of Fig. |4] which refers to the ideal-fluid EoS, 
the comparison of the orbital frequency. The agreement of 
the results of the two codes is excellent, if one ignores the 
initial spurious signal (related to the spurious gravitational- 
wave content of the initial data, which is rapidly propagated 
away). The orbital frequency il is computed in postprocessing 
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The biggest real (i.e. not related to the noise) difference be- 
tween the two curves is during the merger, at around 7.5 ms 
and it is of about 10%, which is consistent with the results 

of a. 

In the right panel of Fig. |4] we plot the amplitude of waves 
as a function of the frequency. The inset shows that the error 
on the amplitude is always at most 10%, which is of the same 
order of magnitude of the one found in the comparison of nu- 
merical codes in binary-black-hole simulations 101 and pro- 
vides here an important consistency check on the numerical 
accuracy and validity of the waveforms of both Whisky and 
SACRA. The discussion of whether, as in [7], also for binary- 
NS-merger waveforms this discrepancy is relevant or not for 
data analysis (namely whether current detectors can or cannot 
distinguish between the waveforms of the two codes) is left to 
a future work, now in preparation II18I1 . 

Finally, in order to give also a strong visual support to the 
goodness of the consistence of gravitational waves computed 
from the two codes, in Fig. |5] which refers to the piecewise- 
polytropic EoS, we show the (/i+)22 waveforms, together 
with the curve predicted by the Taylor- T4 post-Newtonian ap- 
proximation lll8ll67ll68ll . These are the raw data, in the sense 
that no phase shift is performed to achieve the best alignment 
of gravitational waves. The latter procedure is often success- 
fully performed in data-analysis related work and will be in- 
cluded in our future work II18I1 . Nevertheless the similarity of 
the numerical waveforms (both between themselves and with 
respect to the post-Newtonian prediction) in the inspiral part 
is astonishingly good at the highest resolutions adopted here 
(see upper panel of Fig. |5]l. On the contrary, the lower resolu- 
tions (lower panel) are clearly not good enough. 

The situation is somewhat different after the merger. The 
ringdown part shows agreement between the two codes, but 
the waves show some differences both in amplitude and fre- 
quency in the interval after the merger and before the ring- 
down. This is due to the differences in the EoSs, as explained 
in Sec. Ill D 21 Namely, in the present simulations SACRA 
added a thermal part to the piecewise-polytropic EoS, while 
Whisky did not. As shown by the figures, the difference in 
the EoS are irrelevant to the inspiral phase, but not so after the 
merger, as expected. 



IV. CONCLUSIONS 

In this work we have presented the first, detailed com- 
parison of two general-relativistic hydrodynamics codes, the 
Whisky code and the SACRA code. 

We have compared numerical-relativity waveforms and 
other quantities for the last orbits, merger, and collapse of 
equal-mass irrotational binary NS systems, as produced by 
the two independent computer codes. We focused on two ana- 
lytic EoSs, namely the simple ideal-fluid EoS and a piecewise- 
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FIG. 5: (Colour online) Comparison of the waveform (/i+)22. These 
data refer to the piecewise-polytropic EoS. The upper panel refers to 
the "higher resolution" and the lower panel to the "lower resolution" 
(see text for details). As in Fig. [4] the curve labeled T4 is the Taylor- 
T4 post-Newtonian approximation Ila . l67ll6al . 



polytropic EoS, for which we additionally presented more res- 
olutions. The purpose was to perform a stringent consistency 
check of the results from these codes. We found that the wave- 
form frequency and amplitude computed with the two codes 
are in agreement with a discrepancy of at most 10% (this es- 



timate refers to the merger time; the discrepancy is much less 
during the inspiral), which is comparable to the intrinsic error 
of each individual code at the adopted resolutions. We stress 
the fact that this estimated error should be considered here 
an upper limit and that the discrepancy between the waves 
computed in the two codes will be smaller when we will con- 
sider an optimised overlap of the waveforms, in our future 
worklfisll. 

The comparison of purely hydrodynamical quantities, like 
the rest-mass density, shows better results, with a difference 
between the two codes of at most about 1%. This number 
refers however only to global quantities (like maxima and 
norms), but not to point-to-point comparisons, mainly because 
of the small phase difference in the evolution, which makes 
pointwise comparisons meaningless. In fact, even after com- 
pensating for the phase difference, errors larger than 1% are 
seen at some points, noticeably those near the surface of the 
stars. Such errors are related to different implementations of, 
e.g., the atmosphere treatment and do not influence the global 
dynamics in a noticeable way. 

Finally, by comparing other time-dependent spacetime and 
matter quantities, we showed that both codes conserve at high 
accuracy rest mass, energy, and angular momentum, when 
taking into account the emission of gravitational waves. The 
small differences that are present have been related to details 
in the different implementations and grid setups. 

In conclusion, encouraging results have been shown and 
more work is now necessary to assess how the remaining dif- 
ferences in the results may affect the construction of templates 
for gravitational-wave data analysis. This will be the subject 
of a future work |18ll, which may include also more codes in 
the comparison. 
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